function [ result_step_matrix ] = termStep( handles, trace_data, time, voltage)
%termStep is a funciton that calculates the terminal step and returns the
%corresponding step matrix. This step matrix may then be used to ignore the
%terminal step in our other files.

term_reached = [];
term_amp = [];
term_dur = [];
term_iqr =[];
wv_amp = [];
wv_dur = [];
wv_iqr = [];
pre_term_amp = [];
pre_term_dur = [];
pre_term_iqr = [];
processedSteps = [];
term_time = 0;
multi_steps_amp = {};
multi_steps_dur = {};


alpha = handles.alpha;
%filt = exp_filt(trace_data,1,length(trace_data),alpha);
filt = trace_data;
filt(:,2) = time;

% Get the parameters of the step checking from waveformviewer
minStepLength = 0.2*1e-3;
minStepSize = 3;
number_of_steps = 1;
eject_working = 0;


% Get some general stats on the event
wv_amp = mean(trace_data);
wv_iqr = iqr(trace_data);
n = length(filt);

% We will only look for steps that have a minimum length, meaning we don't
% have to start looking for steps until that minimum step length has pass
% and we can stop 'early'

start_point = round(minStepLength/handles.cur_time_tick);
end_point = n - start_point -1;

% This will be a matrix of the data points where there are steps, we will
% be looking in between each of the plautes generated by the steps for
% multiple steps, so to beginning we have a 'step' at the very first point
% and one at the end
%
% The first column will be which data point the step is at the second
% column is the chisquared value of the step and the 3rd column is the
% length ranked value of the step
step_matrix = [1 inf inf; n inf inf];

Chisq=zeros(length(filt)-3,1); %vector of zeros of the length of the data this will be the fit of each of the platueas

% Decide if we should analyze the event for steps or not
% ignore_eject = get(handles.chkbox_ignore_eject, 'Value');
% max_event_length = str2num(get(handles.edt_max_event_length, 'String'));
for i=1:number_of_steps
    
    %Will store the results for the optimal steps within each of the
    %plateaus, since the step matrix contains the start and end points
    %for each of the platues it will have the number of plateaus + 1
    %points, so to store the optimal steps for each of the plateaus we
    %will only need size(step_matrix) - 1
    mini_step_matrix = zeros(size(step_matrix,1)-1,3);
%     disp(['Step Num: ' num2str(i)])
    for j=1:size(step_matrix,1)-1
%         disp(['Mini Plat Num: ' num2str(j)])
        plat_data = filt(step_matrix(j,1):step_matrix(j+1,1),1);
        
        n = length(plat_data);
        
        %start point is always the number of points a minimum step size
        %is
        end_point = n - start_point -1;
        
        if ( (start_point + 1) < end_point)
            
            %where we will hold the chisq data for this plateau
            Chisq=zeros(length(plat_data)-3,1);
            
            %If we have a really big plateau, a speedup
            if n > 2000
                break_size = round(2e-3/handles.cur_time_tick);
                test_points = n:-break_size:1;
                test_points = test_points(end:-1:1);
                left = sum(plat_data(1:test_points(1),1));
                right = sum(plat_data((test_points(1)+1):end,1));
                Chisq_alt = zeros(length(test_points)-3,1);
                
                %     disp('Spacing test_points')
                %     disp(test_points(2)-test_points(1))
                for t=2:length(test_points)-1
                    left = left + sum(plat_data((test_points(t-1)+1):test_points(t),1));
                    
                    dcleft = left/test_points(t);
                    
                    right = right - sum(plat_data((test_points(t-1)+1):test_points(t),1));
                    
                    dcright = right/(n-test_points(t));
                    
                    Chisq_alt(t-1)=(sum((plat_data(1:test_points(t),1)-dcleft).^2)+sum((plat_data((test_points(t)+1):end,1)-dcright).^2))/(n);
                    
                end
                
                [g,best_break]=min(Chisq_alt);
                
                %     disp('first break')
                %     disp(['min= ' num2str(g) ', index= ' num2str(best_break)])
                count = 1;
                while best_break ~= 0 %&& count < 4
                    if best_break+2 > length(test_points) || best_break-2 < 1
                        if best_break-2 < length(test_points)
                            narrow_range = 1:(test_points(best_break+2)-2);
                        else
                            narrow_range = test_points(best_break-2):(length(plat_data)-2);
                        end
                    else
                        narrow_range = test_points(best_break-2):(test_points(best_break+2)-2);
                    end
                    
                    left = sum(plat_data(1:narrow_range(1),1));
                    right = sum(plat_data(narrow_range(1)+1:end,1));
                    %     Chisq_alt = zeros(length(narrow_range)-3,1);
                    
                    narrow_range(1) = [];
                    
                    for t=narrow_range
                        left = left + plat_data(t,1);
                        dcleft = left/t;
                        
                        right = right - plat_data(t,1);
                        dcright = right/(n-t);
                        
                        Chisq(t-1)=(sum((plat_data(1:t,1)-dcleft).^2)+sum((plat_data(t+1:end,1)-dcright).^2))/(n);
                    end
                    
                    %                     disp('begin/end narrow_range')
                    %                     disp([num2str(narrow_range(1)) '/' num2str(narrow_range(end))])
                    ChiInds = 1:length(plat_data)-3;
                    ChiInds(narrow_range-1) = [];
                    Chisq(ChiInds) = inf;
                    [g,h] = min(Chisq);
                    
                    %             disp(['h+1= ' num2str(h+1) ', min= ' num2str(g)])
                    
                    if h+1 == narrow_range(end) || h+1 == narrow_range(1)
                        if h+1 == narrow_range(end) && h+1 ~= length(plat_data)-2
                            best_break = best_break+1;
                        elseif h+1 ~= 2 && h+1 ~= length(plat_data)-2
                            best_break = best_break-1;
                        else
                            best_break = 0;
                        end
                    else
                        best_break = 0;
                    end
                    
                    count = count+1;
                end
            else
                
                % disp(test_points(1:5))
                % disp(test_points(end-5:end))
                left = plat_data(1,1);
                right= sum(plat_data(2:end,1));
                
                %We'll be fancy and do a moving change to the left and
                %right of the step, and then get the dc for each and
                %calculate the chisquared, this is quite a speed up
                %since we only have to do 2 operations as opposed to n
                %operations each time
                
                for t=2:length(plat_data)-2
                    left = left + plat_data(t,1);
                    
                    dcleft = left/t;
                    
                    right = right - plat_data(t,1);
                    
                    dcright = right/(n-t);
                    
                    Chisq(t-1)=(sum((plat_data(1:t,1)-dcleft).^2)+sum((plat_data(t+1:end,1)-dcright).^2))/(n);
                    
                end
            end
            %keyboard
            
            %adjusting the beginning and end points?
            %Chisq(1)=Chisq(2)+1;
            
            %Chisq(length(filt)-1)=Chisq(length(filt)-2)+1;
            
            ChiInds = 1:length(plat_data)-2;
            ChiInds((start_point:end_point)) = [];
            Chisq(ChiInds) = inf;
            
            %     disp(Chisq)
            %         disp(Chisq)
            
            %get the least Chisq
            [g,h]=min(Chisq);
            
            h = h+1;  %Since the first chisq pt represents n=2;
            
            %     disp('h+1')
            %     disp(h)
            % disp('term_point')
            % disp(h)
            
            %figure out the step size
            l=plat_data(1:h,1);
            r=plat_data(h+1:end,1);
            
            leftStepAmp = mean(l);
            rightStepAmp = mean(r);
            leftStepLength = time(h);
            rightStepLength = time(end)-time(h);
            
            stpSize=leftStepAmp-rightStepAmp;
            
            %                 if( (abs(stpSize) < minStepSize) || ( leftStepLength < minStepLength) || (rightStepLength < minStepLength) )
            
            % see if the step meets our critera
            if( ( leftStepLength < minStepLength) || (rightStepLength < minStepLength) || (stpSize < minStepSize) )
                mini_step_matrix(j,:) = [1 inf 0];
                stpSize = 0;
                leftStepLength = 0;
                rightStepLength = 0;
                leftStepAmp = wv_amp;
                rightStepAmp = wv_amp;
                h=1;
            else
                % calculate a rank for the fit based on the pleatu
                % length that we are currently looking at
                rank1=abs(stpSize)*sqrt(length(plat_data));        %expected rel.step accuracy relative to background noise
                rank2=abs(stpSize)*sqrt(length(plat_data)/g);       %measured rel. step accuracy
                h = h + (step_matrix(j,1)-1);               %need to correct h
%                 disp('Step Found')
%                 disp(['point: ' num2str(h)])
%                 disp(['time: ' num2str(time(h)*1000)])
%                 disp(time(h)*1000)
                mini_step_matrix(j,:) = [h g rank1];
                spl=[sqrt(g),h,stpSize,rank1, rank2];      %minimum Chi, index, stepsize,  step*srqt(N), step*sqrt(N/Var)
            end
        else
            mini_step_matrix(j,:) = [1 inf 0];
            stpSize = 0;
            leftStepLength = 0;
            rightStepLength = 0;
            leftStepAmp = wv_amp;
            rightStepAmp = wv_amp;
            h=1;
        end
    end
%     disp('mini_step_matrix')
%     disp(mini_step_matrix)
%     disp('-----------------')
    [maxsig, newstep] = max(mini_step_matrix(:,3));
    
    if( maxsig == 0)
        break;
    else
        step_matrix = [step_matrix; mini_step_matrix(newstep,:)];
        step_matrix = sortrows(step_matrix, 1);
%         disp('step_matrix')
%         disp(step_matrix)
%         disp('------------')
    end
end
%     set(handles.sweep_plot,'YLim',yaxis);
%     ydata = [fin_avg fin_avg];
%     xdata = get(handles.sweep_plot,'XLim');
%     line(xdata,ydata,'color','g');

result_step_matrix = step_matrix;


end

